#include "decomposition.h"
using namespace cv;
using namespace std;

/**
 * decomposition 
 *			get the m, h, d layer
 * @param 
 * @return
 */
void decompose(Mat src, Mat mask, Mat red, Mat delta, Mat dest, double argv2, double argv3, double lambda, double lambda2, int step)
{
	Mat srcback=Mat(cvSize(src.cols, src.rows), CV_8UC3);
	src.copyTo(srcback);

	int arraylength = src.rows * src.cols;

	double (*colordata)[3] = new double[arraylength][3];
	double (*diff)[3] = new double[arraylength][3];
	double (*dered)[3] = new double[arraylength][3];
	double** pcaData;
	long length=src.cols*src.rows;
	double avec[3]={0,0,0};

	// calculate log space
	int i,j;
	long s=0,t=0; 
	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			double r0=log(src.at<Vec3b>(i,j)[0]+1);
			double r1=log(src.at<Vec3b>(i,j)[1]+1);
			double r2=log(src.at<Vec3b>(i,j)[2]+1);	
			colordata[s][0]=r0;
			colordata[s][1]=r1;
			colordata[s][2]=r2;
			avec[0]+=r0;
			avec[1]+=r1;
			avec[2]+=r2;
			s++;
		}
	}

	avec[0]=avec[0]/(double)s;
	avec[1]=avec[1]/(double)s;
	avec[2]=avec[2]/(double)s;

	//subtract the mean vector
	Mat input(length,3,CV_32FC1);
	for(int i=0; i<length; i++)
	{
		for(int j=0; j<3; j++)
			input.at<float>(i,j)=colordata[i][j]-avec[j];
	}

	// PCA and dump the results
	Mat output(length,2,CV_32FC1);

	int numPrincipalComponents = 2;
	PCA pca(input, Mat(), CV_PCA_DATA_AS_ROW, numPrincipalComponents);
	pca.project(input, output);

	int headcount=0;

	for(i=0.15*src.cols;i<0.85*src.cols;i++)
		for(j=0.20*src.rows;j<0.90*src.rows;j++)
			if((i*src.cols+j)%step==0&&mask.at<uchar>(i,j)!=0)
				headcount+=1;

	pcaData = mat_create(headcount, 2);
	int k = 0;
	for(i=0.15*src.cols;i<0.85*src.cols;i++){
		for(j=0.20*src.rows;j<0.90*src.rows;j++)
			if((i*src.cols+j)%step==0&&mask.at<uchar>(i,j)!=0)
			{
				pcaData[k][0] = output.at<float>(i*src.cols+j,0);
				pcaData[k][1] = output.at<float>(i*src.cols+j,1);
				k++;
			}
	}

	printf("%lf, %lf, %lf, %lf, %d\n", pcaData[0][0], pcaData[0][1], pcaData[1][0], pcaData[1][1], headcount);

	// ICA
	double **K, **W, **A, **S;
	int rows, cols, compc;
	compc = 2;
	rows = headcount;
	cols = 2;
	W = mat_create(compc, compc);
	A = mat_create(compc, compc);
	K = mat_create(cols, compc);
	S = mat_create(rows, cols);
	fastICA(pcaData, rows, cols, compc, K, W, A, S);

	double A0,B0,C0,D0;
	A0 = -A[0][0];
	B0 = -A[1][0];
	C0 = -A[0][1];
	D0 = -A[1][1];

	printf(" %lf, %lf\n%lf, %lf\n%d, %d\n", A0, B0, C0, D0, rows, cols);

	Mat resynA(length,1,CV_32FC1);
	Mat resynB(length,1,CV_32FC1);

	double minA=4096,minB=4096;

	for(i=0;i<length;i++)
	{
		double p=output.at<float>(i,0),q=output.at<float>(i,1);     
		resynB.at<float>(i)=(p/A0-q/C0)/(B0/A0-D0/C0);
		resynA.at<float>(i)=(p-B0*resynB.at<float>(i))/A0;
		if(resynB.at<float>(i)<minB)
		{
			minB=resynB.at<float>(i);
		}
		if(resynB.at<float>(i)<minA)
			minA=resynA.at<float>(i);
	}
	pca.backProject(output, input);


	//extract the third layer
	double cmin=255,cmax=0;
	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				//
				double r=input.at<float>(i*src.cols+j,t);
				double s=srcback.at<Vec3b>(i,j)[t];
				colordata[i*src.cols+j][t]=-(exp(r+avec[t])-1)+s;
				if(colordata[i*src.cols+j][t]>cmax)
					cmax=colordata[i*src.cols+j][t];
				if(colordata[i*src.cols+j][t]<cmin)
					cmin=colordata[i*src.cols+j][t];  
				dered[i*src.cols+j][t]=(exp(r+avec[t])-1);
			}
		}
	}

	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				red.at<Vec3b>(i,j)[t]=(colordata[i*src.cols+j][t]-cmin)/(cmax-cmin)*256;
			}
		}
	}

	Mat output2(length,2,CV_32FC1);
	output.copyTo(output2);

	for(i=0;i<length;i++)
	{
		double p=resynA.at<float>(i),q=resynB.at<float>(i);
		double b=0;
		p=argv2*p;
		q=argv3*q+b;
		output.at<float>(i,0)=A0*p+B0*q;
		output.at<float>(i,1)=C0*p+D0*q;
	}


	pca.backProject(output, input);

	double max=0,min=255;
	for(i=0;i<length;i++)
	{
		for(j=0;j<3;j++)
		{
			if(input.at<float>(i,j)>max)
				max=input.at<float>(i,j);
			if(input.at<float>(i,j)<min)
				min=input.at<float>(i,j);
		}
	}

	double rmin=255,rmax=0;
	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=input.at<float>(i*src.cols+j,t);

				r=exp(r+avec[t])-1; //-> the third layer
				if(r<rmin)
					rmin=r;
				if(r>rmax&&r<400)
					rmax=r;
				colordata[i*src.cols+j][t]=r;
			}
		}
	}

	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=(colordata[i*src.cols+j][t]-rmin)/(rmax-rmin)*256;
				dest.at<Vec3b>(i,j)[t]=r;
			}
		}
	}


	double dmax=-1000,dmin=1000;
	for (i=0;i<dest.rows;i++)
	{
		for(j=0;j<dest.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=-colordata[i*dest.cols+j][t]+dered[i*dest.cols+j][t];
				diff[i*dest.cols+j][t]=r;
				if(r<dmin)
					dmin=r;
				if(r>dmax)
					dmax=r;
			}
		}
	}
	if(dmin<-DD) dmin=-DD;
	if(dmax>DD) dmax=DD;

	for (i=0;i<dest.rows;i++)
	{
		for(j=0;j<dest.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=diff[i*dest.cols+j][t];
				delta.at<Vec3b>(i,j)[t]=(r-dmin)/(dmax-dmin)*256;
			}
		}
	}

	// delayer
	printf("start delayer\n");
	Mat gamma = Mat(cvSize(src.cols, src.rows), CV_8UC3);
	Mat gamma2 = Mat(cvSize(src.cols, src.rows), CV_8UC3);
	delayer(delta, gamma, DELAYER_LOW, 15, 0.2);
	delayer(delta, gamma2, DELAYER_HIGH, 15, 0.2);
	printf("end delayer\n");

	for(i=0;i<length;i++)
	{
		double p=resynA.at<float>(i),q=resynB.at<float>(i);
		double b=lambda*gamma.at<Vec3b>(i/src.cols,i%src.cols)[0]+lambda2*gamma2.at<Vec3b>(i/src.cols,i%src.cols)[0];
		p=argv2*p;
		q=argv3*q+b;
		output2.at<float>(i,0)=A0*p+B0*q;
		output2.at<float>(i,1)=C0*p+D0*q;
	}

	pca.backProject(output2, input);

	max=0,min=255;
	for(i=0;i<length;i++)
	{
		for(j=0;j<3;j++)
		{
			if(input.at<float>(i,j)>max)
				max=input.at<float>(i,j);
			if(input.at<float>(i,j)<min)
				min=input.at<float>(i,j);
		}
	}

	rmin=255,rmax=0;
	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=input.at<float>(i*src.cols+j,t);

				r=exp(r+avec[t])-1; //-> the third layer
				if(r<rmin)
					rmin=r;
				if(r>rmax&&r<400)
					rmax=r;
				colordata[i*src.cols+j][t]=r;
			}
		}
	}

	for (i=0;i<src.rows;i++)
	{
		for(j=0;j<src.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=(colordata[i*src.cols+j][t]-rmin)/(rmax-rmin)*256;
				dest.at<Vec3b>(i,j)[t]=r;
			}
		}
	}

	dmax=-1000,dmin=1000;
	for (i=0;i<dest.rows;i++)
	{
		for(j=0;j<dest.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=-colordata[i*dest.cols+j][t]+dered[i*dest.cols+j][t];
				diff[i*dest.cols+j][t]=r;
				if(r<dmin)
					dmin=r;
				if(r>dmax)
					dmax=r;
			}
		}
	}
	if(dmin<-DD) dmin=-DD;
	if(dmax>DD) dmax=DD;

	for (i=0;i<dest.rows;i++)
	{
		for(j=0;j<dest.cols;j++)
		{
			for(t=0;t<3;t++)
			{
				double r=diff[i*dest.cols+j][t];
				delta.at<Vec3b>(i,j)[t]=(r-dmin)/(dmax-dmin)*256;
			}
		}
	}

	delete colordata;
	delete diff;
	delete dered;
	mat_delete(pcaData, rows, cols);
}
